*-------------------------------------------------------------------------------------------------------*
*RESEARCHERS:		Mette Foged, Linea Hasager and Vasil Yasenov
*PROJECT:			Meta-analysis
*DESCRIPTION:		1  Create kernel density plots of effect sizes (full sample, relative and average wages) WEIGHTED
*-------------------------------------------------------------------------------------------------------*


clear all
set more off


use "${data}\database_clean.dta", replace

foreach type in wage empl {
preserve
keep if `type'==1

gen avg=1 if variation=="Area"
replace avg=0 if variation=="Skill" | variation=="Mixture"


*Top and bottom code
local bound 3
gen estimate_coded = estimate
replace estimate_coded = -`bound' if estimate < -`bound'
replace estimate_coded = `bound' if estimate > `bound'
label var estimate_coded "Estimate"

*weights (must be integer, so I multiply by xx)
bysort paper_id: gen p_N=int(1/_N*1000000000)


*remember to generate weights if calculating statistics
qui sum estimate, d
local MEAN_=round(r(mean), .01)
local MED_=round(r(p50), .01)
local SD_=round(r(sd), .01)
local N_=r(N)

qui sum estimate if avg==1, d 
local N_avg=r(N)
local MEAN_avg=round(r(mean), .01)
local MED_avg=r(p50)



qui sum estimate if avg==0, d
local N_rel=r(N)
local MEAN_rel=round(r(mean), .01)
local MED_rel=r(p50)



*redefine weights (must be integer, so multiply )
bysort paper_id avg: gen p_N_a=int(1/_N*1000000000)
gen p_score=int(1/Score*1000000000)
gen p_se=int(1/se*1000000000)



local bound 3
*Average vs. relative wages	, estimates weighted equally
twoway (kdensity estimate_coded if avg==0, width(0.2) xline(0, lcolor(black)) color(none) lcolor(black) lpattern(longdash)  lwidth(medthick) xlabel(-`bound'(.5)`bound') ylabel(0(0.5)1.5)) (kdensity estimate_coded if avg==1 , width(0.2) color(gs12%30) lcolor(black) lpattern(solid) xlabel(-`bound'(.5)`bound') ylabel(0(.5)1.5) xtitle(`"{fontface Times New Roman:Estimate}"') ytitle(`"{fontface Times New Roman:Density}"')) ///
, graphregion(color(white)) legend(off) saving(estimates, replace) 
graph export  "${out}\histogram_avg_rel_`type'.eps", replace 



local bound 3
*Average vs. relative wages	, papers weighted equally
qui sum estimate [weight=p_N_a] if avg==1, d
local MEAN_avg=round(r(mean), .01)
local MED_avg=round(r(p50), .01)
local N_avg=r(N)
qui sum estimate [weight=p_N_a] if avg==0, d
local MEAN_rel=round(r(mean), .01)
local MED_rel=r(p50)
local N_rel=r(N)

twoway (kdensity estimate_coded if avg==0 [fweight=p_N_a], width(0.2) xline(0, lcolor(black)) color(none) lcolor(black) lpattern(longdash)  lwidth(medthick)  xlabel(-`bound'(.5)`bound') ylabel(0(.5)1.5)) ///
(kdensity estimate_coded if avg==1 [fweight=p_N_a], width(0.2) color(gs12%30) lcolor(black) lpattern(solid) xlabel(-`bound'(.5)`bound') xtitle(`"{fontface Times New Roman:Estimate}"') ytitle(`"{fontface Times New Roman:Density}"')) , graphregion(color(white)) legend(off) saving(papers, replace)  
graph export  "${out}\histogram_avg_rel_paper_`type'.eps", replace

local bound 3
*Average vs. relative wages	, weighted by score
qui sum estimate [weight=p_score] if avg==1, d
local MEAN_avg=round(r(mean), .01)
local MED_avg=round(r(p50), .01)
local N_avg=r(N)
qui sum estimate [weight=p_score] if avg==0, d
local MEAN_rel=round(r(mean), .01)
local MED_rel=round(r(p50), .01)
local N_rel=r(N)

twoway (kdensity estimate_coded if avg==0 [fweight=p_score], width(0.2) xline(0, lcolor(black)) color(none) lcolor(black) lpattern(longdash)  lwidth(medthick)  xlabel(-`bound'(.5)`bound') ylabel(0(.5)1.5)) (kdensity estimate_coded if avg==1 [fweight=p_score], width(0.2) color(gs12%30) lcolor(black) lpattern(solid) xlabel(-`bound'(.5)`bound') xtitle(`"{fontface Times New Roman:Estimate}"') ytitle(`"{fontface Times New Roman:Density}"')) , graphregion(color(white)) legend(off) saving(score, replace) 
graph export  "${out}\histogram_avg_rel_score_`type'.eps", replace

local bound 3
*Average vs. relative wages	, weighted by SE
qui sum estimate [weight=p_se] if avg==1, d
local MEAN_avg=round(r(mean), .01)
local MED_avg=round(r(p50), .01)
local N_avg=r(N)
qui sum estimate [weight=p_se] if avg==0, d
local MEAN_rel=round(r(mean), .01)
local MED_rel=round(r(p50), .01)
local N_rel=r(N)

twoway (kdensity estimate_coded if avg==0 [fweight=p_se], width(0.2) xline(0, lcolor(black)) color(none) lcolor(black) lpattern(longdash)  lwidth(medthick)  xlabel(-`bound'(.5)`bound') ylabel(0(.5)1.5)) (kdensity estimate_coded if avg==1 [fweight=p_se], width(0.2) color(gs12%30) lcolor(black) lpattern(solid) xlabel(-`bound'(.5)`bound')  xtitle(`"{fontface Times New Roman:Estimate}"') ytitle(`"{fontface Times New Roman:Density}"')) , graphregion(color(white)) legend(off)  saving(se, replace) 
graph export  "${out}\histogram_avg_rel_se_`type'.eps", replace



erase estimates.gph
erase papers.gph
erase score.gph
erase se.gph

restore
}


